library(tidyverse)
library(plotly)
library(sf)
library(tigris)
library(leaflet)
library(censusapi)
library(zoo)
years <- 2017:2020
quarters <- 1:4
types <- c("Electric","Gas")
pge_17_to_20_elec_and_gas <- NULL

for (year in years) {
  if (year == 2020) {
    quarters <- 1:2
  }
    for (quarter in quarters) {
      for (type in types) {
        filename <-
          paste0(
            "PGE_",
            year,
            "_Q",
            quarter,
            "_",
            type,
            "UsageByZip.csv"
        )
        temp <- read_csv(filename)
        pge_17_to_20_elec_and_gas <- bind_rows(pge_17_to_20_elec_and_gas, temp)
        saveRDS(pge_17_to_20_elec_and_gas, "pge_17_to_20_elec_and_gas.rds")
    }
  }
}
ca_counties <- counties("CA", cb=T, progress_bar = F)

bay_county_names <-
  c(
    "Alameda",
    "Contra Costa",
    "Marin",
    "Napa",
    "San Fransisco",
    "San Mateo",
    "Santa Clara",
    "Solano",
    "Sonoma"
  )

bay_counties <-
  ca_counties %>%
  filter(NAME %in% bay_county_names)

usa_zips <-
  zctas(cb = T, progress_bar = F)

bay_zips <-
  usa_zips %>%
  st_centroid() %>%
  .[bay_counties, ] %>%
  st_set_geometry(NULL) %>%
  left_join(usa_zips %>% select(GEOID10)) %>%
  st_as_sf()
pge_final <-
  pge_17_to_20_elec_and_gas %>%
  filter(CUSTOMERCLASS %in% 
           c(
             "Elec- Residential",
             "Elec- Commercial",
             "Gas- Residential",
             "Gas- Commercial"
           ),
        ZIPCODE %in% bay_zips$ZCTA5CE10) %>%
  mutate(
    TOTALKBTU = case_when(
      !is.na(TOTALKWH) ~ TOTALKWH * 3.412,
      !is.na(TOTALTHM) ~ TOTALTHM * 99.976
    )
  ) %>%
  group_by(YEAR, MONTH, CUSTOMERCLASS) %>%
  summarize(
    MONTHLYKBTU = 
      sum(
        TOTALKBTU,
        na.rm = T
      )
  ) %>%
 mutate(Date = as.yearmon(paste(YEAR, MONTH, sep = "-")))
pge_chart <-
  pge_final %>%
  ggplot() +
  geom_line(
    mapping = aes(
      x= Date, 
      y= MONTHLYKBTU, 
      color = CUSTOMERCLASS), 
      stat = "identity" 
      ) +
  labs(
    y = "Monthly kBTU",
    title = "PG&E Bay Area Monthly Electricity and Gas Usage, 2017-2020",
    fill = "Electricity Type"
  ) 

  pge_chart %>% ggplotly()

This line chart shows observable changes in energy consumption that may be attributable to the COVID-19 pandemic. I compared electricity and gas usage from April-June of 2019 to April-June of 2020 and found that for electricity usage, residential consumption is higher in April-June 2020 than in April-June 2019. Commercial consumption follows the opposite trend — it’s higher in April-June 2019 than in April-June 2020. The largest difference in commercial electricity consumption happens in April, representing a 21.6% decrease from 2019 to 2020, while the largest difference in residential electricity consumption happens in May, representing a 15% increase from 2019 to 2020. This suggests that the first full month of California’s stay-at-home-order (April) had a greater affect on businesses and other commercial buildings, whereas the change in residential electricity consumption peaked in May.

Overall, gas consumption followed a similar trend as electricity consumption, but there are a couple of exceptions. In April of 2020 commercial gas consumption was actually 3.3% higher than in April of 2019, and in May of 2020 residential gas use was 14% lower than in May of 2019. The largest difference in residential gas use occurred in April, when consumption was 31% higher in 2020 than in 2019. Moreover, the largest difference in commercial gas use occurred in May, when consumption decreased ~17% from 2019 to 2020. This trend suggests the opposite conclusion from that of electricity consumption; the first full month of the stay-at-home order had a significant impact on residential gas consumption, whereas commercial gas consumption experienced its greatest change in May.

We cannot definitively attribute these trends to COVID-19, but it’s evident that there are changes in consumption and it’s reasonable to assume that the pandemic is, to a certain degree, playing a role in causing these differences.

Lastly, it’s important to recognize the differences in consumption data for September 2017 from the rest of the trend. It appears that there’s an error in the PG&E data, which is leading to the appearance of spikes in energy consumption for that month. If the data were correct, we could expect to see a trend similar to that of the other years.

pge_res_elec_covid <-
  pge_17_to_20_elec_and_gas %>%
      filter(YEAR %in% (2019:2020), MONTH %in% (4:6), 
      CUSTOMERCLASS == "Elec- Residential") %>%
      mutate(
        ZIPCODE = ZIPCODE %>% as.character()
      ) %>%
      group_by(ZIPCODE, YEAR) %>%
      summarize(
       TOTALKWH = 
         sum(TOTALKWH, na.rm = T)) %>%
       pivot_wider(
         names_from = YEAR,
         values_from = TOTALKWH
       ) %>%
     rename(
       KWH2019 = "2019", KWH2020 = "2020"
     ) %>%
     mutate(
     percent_change = 
       ((KWH2020 - KWH2019)/KWH2019) * 100
     ) %>%
     right_join(
     bay_zips %>% select(GEOID10),
     by = c("ZIPCODE" = "GEOID10")
     ) %>% 
     filter(KWH2019 != 0, KWH2020 !=0, !is.na(KWH2019), !is.na(KWH2020)) %>%
     st_as_sf() %>% 
     st_transform(4326)
pge_comm_elec_covid <-
  pge_17_to_20_elec_and_gas %>%
      filter(YEAR %in% (2019:2020), MONTH %in% (4:6), 
      CUSTOMERCLASS == "Elec- Commercial") %>%
      mutate(
        ZIPCODE = ZIPCODE %>% as.character()
      ) %>%
     group_by(ZIPCODE, YEAR) %>%
     summarize(
       TOTALKWH = 
        sum(TOTALKWH, na.rm = T)
     ) %>%
    pivot_wider(
       names_from = YEAR,
       values_from = TOTALKWH
     ) %>%
    rename(
      KWH2019 = "2019", KWH2020 = "2020"
    ) %>%
    mutate(
    percent_change = 
      ((KWH2020 - KWH2019)/KWH2019) * 100
    ) %>%
    right_join(
    bay_zips %>% select(GEOID10),
    by = c("ZIPCODE" = "GEOID10")
    ) %>% 
    filter(KWH2019 != 0, KWH2020 !=0, !is.na(KWH2019), !is.na(KWH2020)) %>%
    st_as_sf() %>% 
    st_transform(4326)
res_pal <- colorBin(
  palette = "PiYG",
  domain = 
    pge_res_elec_covid$percent_change,
  bins = c(-75, -20, -10, 0, 10, 20, 30)
  )

leaflet() %>% 
  addTiles() %>% 
  addPolygons(
    data = pge_res_elec_covid,
    fillColor = ~res_pal(percent_change),
    color = "white",
    opacity = 1,
    fillOpacity = 1.5,
    weight = 1,
    label = ~paste0(
      round(percent_change), 
      " percent change in ",
      ZIPCODE
    ),
    highlightOptions = highlightOptions(
      weight = 2,
      opacity = 1
    )
  ) %>%
   addLegend(
    data = pge_res_elec_covid,
    pal = res_pal,
    values = ~percent_change,
    title = "Percent Change in <br> April-June Residential <br> kWH,     2019-2020"
  ) 

I chose to represent “before COVID began” as cumulative electricity consumption from April-June of 2019, and “after COVID began” as cumulative electricity consumption from April-June of 2020. I chose to compare 2019 to 2020 because I believe 2019 data would provide the most relevant information on energy use in the Bay Area as a baseline year. I chose the months April-June because April was the first full month of California’s shelter-in-place, and because June is the most recent available month for 2020 data. While the effects vary by neighborhood, with some zipcodes experiencing at least a 20% increase in energy use and others experiencing a decrease in energy use, the majority of zipcodes have experienced an increase in residential electricity use since the pandemic began.

res_pal <- colorBin(
  palette = "RdBu",
  domain = 
    pge_comm_elec_covid$percent_change,
  bins = c(-90, -70, -50, -30, -10, 0, 25, 50, 125, 200, 300)
)

leaflet() %>% 
  addTiles() %>% 
  addPolygons(
    data = pge_comm_elec_covid,
    fillColor = ~res_pal(percent_change),
    color = "white",
    opacity = 1.5,
    fillOpacity = 2,
    weight = 1,
    label = ~paste0(
      round(percent_change), 
      " percent change in ",
      ZIPCODE
    ),
    highlightOptions = highlightOptions(
      weight = 2,
      opacity = 1
    )
  ) %>%
   addLegend(
    data = pge_comm_elec_covid,
    pal = res_pal,
    values = ~percent_change,
    title = "Percent Change in <br> April-June Commercial <br> kWH, 2019-2020"
  ) 

My choice in defining “before COVID began” and “after COVID began” follows the same logic as above. Similar to the percent change in residential energy use, the percent change in commercial energy use also varies by neighborhood. However, as is clear by the map, the overwhelming majority of zipcodes have experienced a decrease in commercial electricity use since COVID began, some even as large as 84%. There are a few zipcodes, shown in dark blue, that indicate a large increase in commercial energy use. While we cannot definitively attribute the cause of these increases, it’s possible that these commercial buildings provide services that are in higher demand during the pandemic, such as medical care. Even with these few exceptions, overall commercial electricity use has decreased across the Bay Area since COVID began.